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Abstract 

This paper describes the design and testing of an 
Extended Kal«an Filter (EKF) for ground attitude 
determination, misalignment estimation and sensor 
calibration of the Earth Radiation Budget Satellite 
{ERBS). Attitude is represented by the quaternion of 
rotation and the attitude estimation error is defined 
as an additive error. Quaternion normalization is used 
for increasing the convergence rate and for minimizing 
the need for filter tuning. The paper presents the 
development of the filter dynamic model, the gyro error 
model and the measurement models of the Sun sensors, 
the IR horizon scanner and the magnetometers which are 
used to qenerate vector measurements. The filter is 
applied to real data transmitted by ERBS se nso ™- 
Results are presented and analyzed and tne tPtr 
advantages as well as sensitivities are discussed. On 
the whole the filter meets the expected synergism, 
accuracy and robustness. 


I. INTROOUCTION 

An important part of spacecraft ground support is 
attitude determination, sensor alignment, and sensor 
calibration. In the past, at Goddard Space Flight 
Center (GSFC) in the Flight Dynamics Division (FDD) 
each task was performed separately, usually using a 
relatively small state . The use of more sophisticated 
algorithms has been suggested in the literature, but 
they have not yet been tested with real spacecraft data 
for ground processing in Flight Dynamics. 

The purpose of this study was to design and test an 
Extended Kalman Filter (EKF). The filter was designed 
for the Earth Radiation Budget Satellite (ERBS). ERBS 
is equipped with the following sensors which are used 
for attitude determination: 2 redundant Inertial 
Reference Units (IRUs) each containing 3 single-axis 
gyroscopes, 2 digital fine Sun sensors (FSSs), 2 Infra- 
red (IR) horizon scanners, and 1 three-axis 
magnetometer. The state estimated by the filter 
consists of the attitude parameters (quaternion), 
sensor misalignments for the Sun sensor, magnetometer 
and gyros, biases for the Sun sensor, horizon scanner, 
magnetometer, and gyros, and scale factor corrections 

for the Sun sensor, magnetometer, and gyros. The 
filter was tested using real spacecraft data 
transmitted to Earth by ERBS. 

Kalman filters have not been used for ground 
attitude processing In the FDD at GSFC. The current 
ground support software implements single frame and 
batch estimators and, as mentioned before, much of the 
calibration effort Is performed separately from the 
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attitude determination. The EKF designed for ERBS 
allows for all of the calibration to be performed along 
with the attitude determination. 

The use of the extended Kalman filter (EKF) for 
spacecraft attitude determination has been dealt with 
quite extensively in the past. Kau et al . [1] as well 
as Farrell [2], for example, used an ad-hoc solution to 
the problem of estimating the Euler angles directly 
from vector measurements. A more general approach to 
this problem was presented in [3]. The problem of 
estimating the direction cosine matrix directly from 
vector measurements was discussed in [4]. The filter 
which was required there was a linear one with some 
adaptation. A general analytic exposition of the use of 
the EKF for spacecraft attitude determination was given 
by Lefferts, Markley and Shuster [5]. Reference [6] 
dealt with the problem of estimating the attitude 
quaternion from vector measurements. Basically, the 
estimated quantity was the difference between the best 
known value of the quaternion and Its true value. This 
difference was defined as a four component add It to 
quantity. Because of this definition the estimate of 
the quaternion is not necessarily normal unless it 
converges to the correct quaternion. It was found that 
normalization of the estimated quaternion during the 
filtering process speeds up convergence and eliminates 
the need for filter tuning. In other references, e.g. 

[5] , [7] and for on board attitude determination 
software which Is used in LANDSAT 4 and is planned to 
be used in the GRO and EP spacecraft a mul tipi icatto 
quaternion difference is used. Since it is assumed that 
this difference quaternion is small and as for small 
rotations the scalar part of the quaternion is close to 
1, those algorithms are estimating only three attitude 
- error components. Obviously, estimation of an 
additive quaternion error of four parameters plus the 
induced normality constraint is equivalent to 
estimating three parameters. Because of our good 
experience with the additive quaternion error approach 

[6] we chose to Implement this approach in the present 
EKF algorithm. 

In the next section we Introduce the algorithms 
developed for the ERBS EKF. 


II. THE EXTENDED KALMAN FILTER ALGORI THM 

The EKF algorithm is based on the following assumed 
models: System model : 



X - f(X(t),t) + *(t) 

(2.1) 

Measurement Model: 


where: 

ik " hk(X(^k)) + ^k 

(2-2) 


X(t) * state vector. 
y( t)« zero mean white process, 
y^- zero mean white sequence. 


The EKF algorithm is as follows [8], The 
update of the state estimate and of the 
error covariance are performed as follows: 

measurement 

estimation 
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where: 


State Estimation Update : 

2 k (+) - $ k (-) ♦ K k (i k - Jj k (l k (-))] (2.3) 

where the Gain Matrix Is evaluated as follows: 

K k ■ p k (-)H k T (i k (-))[H k (i k (-))P k (-)H k T (i k (-)) + R ^* 1 

... ( 2 . 4 ) 

Error Covarlince Update : 


Y - effective measurements. 

M A j • transformation matrix from the nominal (non- 
mlsal Igned) sensor to body coordinates, 
w T\meas " unit vector as measured by the sensor In the 
A sensor misaligned coordinates. 

A(fl) • transformation matrix from the inertial to the 
body coordinates as a function of the 
estimated quaternion. 

- the measured unit vector as known In the 
Inertial coordinates. 


p k {+) - [I - K k H k ]P k {-)[I - K k H k ] T + K k R k K k T (2.5) 

The propagation of the state estimate and the error 
covariance are accomplished using: 

Silt £ Estl ruatlw Propagation : 


i(t) - f(itt).t) 

(2.6) 

Error Covtr lance Propagation: 



p(t) * F(i(t),t)P(t) + P(t)F T (i(t).t) + Q(t) 

(2.7) 

where 

F(X(t),t) 

*i(t) 

1 

1 i(t)-i(t) 

(2.8a) 

■>h( i(t)) 

H(X(-)) - 

^i(t) 


(2.8b) 


While the traditional EKF algorithm updates the state 
estimate according to (2.3), we use y (as computed In 
(2.10)) to update the state estimate as follows 


i k (+) - X k (-) + K k X k (2.11) 


To reconcile this apparent deviation from the ordinary 
EKF algorithm, define d i k as follows 

d2 k - Z k - h k (£ k <-)) (2.12) 

then (2.3), the state update equation in the ordinary 
EKF algorithm, reads 


£k(+) - £*(-) + K k di k (2.13) 

Next define *(t k ) as 

£{t k ) - £ k (-) + *(t k ) (2.U) 

expand (2.2) In a Taylor series expansion about & k (-) 
and omit terms of second and higher order of x(ti,). 
This yields k 


R k • covariance matrix of white sequence. 
Q - spectral density matrix of g(t). 


The EKF rather than the linear KF algorithm must be 
used because the measurement vectors obtained from the 
sensors are non-linear functions of the state vector. 
The state vector was selected to be: 


fl | 

*■1 

*«! 


£i l 

*1 

*■! 

SJ 


4 quaternion components 
3 gyro scale factor errors 
6 gyro misalignment angles 
3 gyro biases 
3 FSS misalignment angles 
2 FSS scale factor errors 
2 FSS biases 

2 IR horizon scanner biases 

3 magnetometer scale factors 
6 magnetometer misalignments 


(2.9) 


i k ■ ll k <i k (-)) + H k S(t k ) + y k (2.15) 

where H k Is as defined In (2.8b). When from (2.15) 
Is substituted Into (2.12) we obtiln 


d ik ■ H kX (t k ) + y k (2.16) 


thet Is, di k JH linearly related la x(t k ). An 
inspection of (2.13) reveals that the EKF estimates 
i(t k ), which according to (2.16) Is linearly related to 
the Affect yve measurement di k , and then adds the 
estimate, A(t k ), to X k (-), the best estimate of 
l(t k ). As will oe seen in the ensuing, also our use of 
y, as defined In (2.10), In the state update equation, 


(2.11), amounts to estimating y(t k ), which 1s A linearly 
related to y . and adding the estimate to Xjc(-). In 
fact, to show the latter we only have to snow that 
£(t k ) Is linearly related to y. This will Indeed be 
shown In Section IV. 


HI. IH£ PYNAMICS MODEL 

The states which vary in time are the attitude 
parameters and bias states which are modeled as Markov 
rather than as bias states. (The reason for this 
modeling will be discussed later). The scale factors 
and misalignments are assumed to be constant in time. 



3 magnetometer biases 


The attitude matrix is given in terms of the 
quaternion, q, as follows 


Following the tradition of the NASA Goddard’s Flight 
Dynamics Division we used vector measurements to update 
the EKF. (It should be noted that this Is not a must 
but rather a choice). The effective measurements which a 
are used to update the filter are defined as follows 

* ■ M ATtfT\me»s * ( 2 - 10 ) 


I qf-qf-qf+qi 

j 2(q!q 2 -q3q 4 ) 
j 2(qjU3+<l2 c l4) 


2 ( QiQ2 +< l3 £ l4 ) 
2 ( < l2 c l3- t ll<l4) 


I 

2{qiq3-q 2 q 4 )j 
2(q2q3 + qiq 4 ) l 
-qf -q|+q|+q| I 


(3.1) 
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The quaternion changes in time according to [8, pp. 
511, 512] 


a ■ 2q 


(3.2) 


where: 


«z " w y "x 


-w z 0 


"y -W x 


(3.3) 


and where w„ v 7 are the components of the spacecraft 
angular velodity vector resolved in the spacecraft 
coordinates. The true quaternion of the spacecraft 
propagates in time according to (3.2). We cannot 
compute q precisely since we do not know precisely the 
initial quaternion nor do we know w precisely as it is 
a measured vector and the measurement contains errors. 
The measured angular velocity can be written as 


where 


S 

a 

dw 


w - y + dw 


gyro reading. 

true angular velocity. 

vector of gyro errors. 


(3.4) 


Since the true quaternion propagates according to (3.2) 
we propagate the estimated quaternion in a similar 
manner; that is, we propagate it according to 


Ofl 


(3.5) 


where ft has the form of (3.3) but its elements are the 
elements of the measured angular rate §. Now a matrix 
dft can be defined as follows 

ft - ft - dft (3.6) 

Substitution of (3.6) Into (3.2) results in 

g * Oft - dftg (3.7) 

When (3.5) is subtracted from (3.7) we obtain 

g - ft » ft(ft - ft) - dftg (3.8) 

As discussed In the introduction, we define an additive 
quaternion error as follows 


dg - a - g 

Then (3.8) can be written as 

dg - Qdg - dftg 

A matrix, B, can be defined as follows 


1" 

-q 4 

<13 

«2 _ | 

\ 1 

-<13 

-q 4 

'll j 

2 j 

q 2 

-qi 

-q 4 I 

I 

. 1l 

<12 

<13-1 


(3.9) 


(3.10) 


(3.11) 


and used In (3.10). However A since g itself is not 
known, we use its estimate, g to compute (3.11). When 
this is done, we can write (3.10) as follows 


| dg - ftdg + §dw j 


(3,12) 


I. 


where § is computed as in (3.11) using g rather than ft. 
Equation (3.12) is the dynamics equation of the 

additive quaternion error. 

Equation (3.12) cannot be used as a dynamics model 
in an EKF since the vector of gyro errors, dg, is not a 
white noise vector. It could be modeled though as a 
linear system excited by a white noise. Consequently 
this linear model can be augmented with the dynamics 
model of (3.12). The augmented model is linear and Is 
driven by a white noise vector hence the model can 
legitimately used by the EKF (8]. To accomplish that we 
use the following standard gyro error model. 


|3w x | 1$qx 
I I I 9 
|dw y |-|0 

! y l I 

|dw z | |0 

L J L 


s gy 0 


gz 


I V 

I 

l«y 

I y 

i?. 


j° ®gxy ®gxz 


+ [ 9 gy* 0 
j®gzx ®gzy 0 


gyz 


IE, 

I 

♦lb, 

I 

lb, 


gz 


j n glx 

+ jUgJy 

! n giz 


(3.13) 


where 


w gxy* 0 gxz» 0 gyx i 0 gyz 


> 0 gzx» 0 gzyl 


Sg ■ t s gx» s gy* s gz^ 

^ - l 0 . 

kg - [bgj(, bgy, bg j ] 

Ugl * t n glx» n gly’ n g 1 z ^ 

and T denotes the transpose, Sn, and b- are as 
explained In (2.9) and n^i 1* a White 3 noise vSctor. We 
can write (3.13) as follows 


(3.14a) 
(3.14b) 
(3.14c) 
(3 . 14d) 


|3w y | |w y 0 0 w v w 7 0 0 0 0 1 

III ^ „ n 

I dw v | - 1 0 w 0 0 0 w w, 0 0 0 

1,1 

jdw 7 j |0 0 w 7 0 0 0 0 w x w y 0 

i_:i l 


0| 


|"gix 


0 1 H ♦ l^qly 
I 1 9 y 


II 

J 


where 


(Sq T . 9g T - fcg) 


j n g 1 z 
. (3.15) 

(3.16) 


Define the following matrices 
0" 


! w x 0 


u - |0 

I 

|0 

L 


Wy 0 


w - 


"Wy w 2 0 0 0 0 

0 0 


w x w 2 0 0 


0 0 


0 0 w x Wy 


(3.17b) 


... (3.17a) 

then (3.15) can be written as 

dw - (U|W|i] 2 * + n^i 

The vectors and g- contained in x* are constants, 
therefore 


(3.18) 
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(3.19b) 


where 


Sg* 


(3.19a) 


Sg 


The gyro bias vector, tg, may actually be time- varying 
so they are more adequately modeled as Markov states as 
follows [8], 


d 

dt 


1 X 
UQ 

j-l/tg 

0 

*'IH 

| n g2xj 

1 b gy 1 ’ 
1 9 1 

1 0 
1 

-1/tg 

0 1 1 bgy 1 * 

ll 9J | 

j Vyj 

Ib gz ! 

L_l 

1 0 

0 

-1/t | |b | 

:il _i 

j n g2z| 


(3.20) 


where t g Is the time constant of the Markov states and 

Hg2 " f n g2x* n g2y* n g2zl (3.21) 

is the white noise vector which drives the Markov 
states. Define the matrix T g as follows 


*s 

£ 

£ 


0 $y» ®szJ 


£ 


tc A - c b] 


t s mx 

[9. 


•'my* ^mz 


mxy> 0 mxz* 0 myx> 9 myz» G mzx’ e mzyJ 

... (3.28) 

The seven sensor states (of the Sun sensor, IR horizon 
scanner and magnetometer) which are listed in (3.26) 
and in (3.28), are augmented with the quaternion error 
and gyro states to form the attitude augmented state 
vector, This vector is that shown in (2.9) when fl is 
replaced by dfl. The differential equation which governs 
the propagation of £ is obtained by combining the 
linear differential equations of the components of the 
attitude augmented state vector. Accordingly the 
augmentation of (3 .n), (3.18), (3.19), (3.23), (3.25) 
and (3.27) yields 


-1/tg 0 


-v t c 


0 

0 

-Vt c 


(3.22) 


then (3.20) can be written as 


fcg ■ T gkg + Bg2 


(3.23) 


The other bias states in the fine Sun sensor, IR 
horizon scanner, and magnetometer which are listed 
In (2.9) and will be mentioned in the development 
of the sensor error models, are also modeled as 
Markov states as follows. Define the following 
matrices 


I -1/t. 


"I 


I “ I 
L o -i/t s j 


(3.24a) 


l-l/th 0 I 
i I 
L o -1/thJ 

... (3.24b) 


: i/t m 0 o' 

0 -i/t m 0 
0 0 -1/t, 


(3.24c) 


then 


d 

dt 


fc s - T s b s + c s 

(3.25a) 

th - T h b h + n h 

(3.25b) 

4n " Vm + 2m 

(3.25c) 


where 

h] - [b A . b B ] 


bS “ ( d r> d pl km * ( b ir 


b m y> b mzl 
... (3.26) 

These vectors denote "biases" as defined in (2.9). The 
scale factor and misalignment states of the sensors 
which also are a part of the state vector listed in 
(2.9), are assumed Constant. That is 




®m 


(3.27) 


21 

< CO 

< QO 

l=r~ 

< CO 

i 

i 

i i i 
! ! ! 

i 

i 

j 

i 

i 

1 

1 

1 

1 

1 '1 

1 dal 
1 | 

1* ~ 
jBflg] 

1 

I i i 

i 

i 


i 

1 


1 

• | 

! i 

.j. 

! ! ! 

i 

i 

i 

1 

1 

* 

1 

j 

1 

.1 


i 

.j. 

! i 

• i * 

i 

i 

i 

i 

1 

h| 

j °gi 

1 

.1 

! ! ! 

i 

.j. 

I I I' 

• i • 

i 

i 

1 

N 

1 

1 

1 

| 

i i I 

i 

.j. 

I i"i 

| 

j 

1 

Yl + 

i f 

1 

1 

1 

• | 

! | j 

i 

.j. 

.MT 

* 1 * 
i 

i 

V 

. I 

1 

h| 


.| 

! ! ! 

i 

.j. 

.LIT?! 

• i • 

i 

i 

1 

1 

1 

h| 

[ n h 

1 

! ! 1 

i 

i i i 

• i * 

i 

V 

’l 

1 H 

l‘"‘ 

"l 

i ! ! 

i 

i'ii 

i 

i 

V 

1 

‘i 

1 

hi 

l"" 

' 1 

i ! I 

V 

V i 'i 

Y 

• 1 * ■ 1 
l T m| 

J ^mj 

.... 
| % 








... (3.29a) 


da | 

• • • I 

.*i 

.M 

.M 

. s ’l 

. £s |- 

.M 

. a ‘! 

5.1 

'«*! 

yi 


which is of the form 


i * F(£)2 + fl (3.29b) 

The spectral density of the elements of the white noise 
driving Markov states in £ is related to the individual 
states they drive according to the well known relation 

iu "t /I 2 / m S T,o where G 1 is the spectral density of 
the white nois§ driving state i, T i is the time 
constant of this Markov state and S i 0 is the initial 

standard deviation of the state. The matrix F($) is the 
one defined in (2.8a). 

The estimation problem dealt with in this paper is 
characterized by a linear dynamics equation. The system 
dynamics Is determined by (3.5), (3.19), (3.23), (3.25) 
and (3.27). It is easy to see that when these equations 
are augmented into one equation we obtain an equation 
of the form 

1 - f (t)X + n (3.30) 
where A is given by (2.9) and f(t) is the following 
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matrix 


f(t) 


I 

I 

I 


(3.31) 


The white noise vector n is of no consequence when 
dealing with the role of (3.30) in the estimation 

process since according to (2.6) the propagation of £ 

requires only the evaluation of f(t). 


IV. IHE MEASUREMENT MQ PEL 

As mentioned in Section II the effecting 
measurements which are used to update the filter are 
defined as follows 

j * ■ ^AT-T’ ,meas * j 

where: 

y - effective measurements. 

- transformation matrix from the nominal (non- 
misaligned) sensor to body coordinates. 

W T , - unit vector as measured by the sensor in the 
’ sensor misaligned coordinates. 

A(fl) » transformation matrix from the Inertial to 
the body coordinates as a function of the 
estimated quaternion. 

y r - the measured unit vector as known In the 

inertial coordinates. 


Consider first W*. The ideal sensor measures in its 
misaligned coordinates the vector Wj » . Since the sensor 
is not Ideal, it adds to the measured vector the error 
term dW-j- * , hence 

-T’ ,meas ’ -T’ + d *T’ (4 ' 4) 
Substitution of (4.4) Into (4.3a) yields 


*A * M AT^T’ + dy T’) 


Now 

h T’A * m T’T h TA 

For small misalignment angles 


where 


From (4.6) 


"T’T 


I + 0 


0 0 Z 
-e z 0 e x 

- e y - 9 x °- 


m AT " m AT’ m T’T 

Substitution of (4.7) into (4.9) yields 


(4.5) 

(4.6) 

(4.7) 

(4.8) 

(4.9) 


H AT - Haj-U + 9) (4.10) 

When (4.10) Is substituted into (4.5) and the term 
containing products of errors is dropped, the following 
is obtained 

Ha “ Mat'Wj. + + M^-j.dW-p (4.11) 

Next we address defined in (4.3b). Using the 
definition of dfl in (3.9) we van write 

A(fl) - A(fl - dfl) (4.12) 

Using Taylor series expansion A(fl) can be approximated 
to within first order terms as follows 

A(fl) - A( fl ) - dq, (4.13) 

i-1 *0i |fl 


substitution of (4.13) into (4.3b) yields 

¥a - ACflIXi dq i (4 ' 14) 


In the ideal (nominal) situation the sensor is well 
aligned and, \n addition, introduces no measurement 
errors. Also, fl, the estimate of a is perfect and is, 
thus, equal to a itself. Therefore, using (4.1), we 
obtain 

1 ■ H A T*T\meas - A (fl)5£l * m AT’Mt- ‘ A (3)Vl * 0 ( 4 -?) 

Any deviation from the nominal will be reflected In y. 
If the deviations are small, then y will be related 
linearly to them. It is our purpose in this section to 
derive the linear relations between the effective 
measurement y and those deviations which are actually 
the error states in & (whose time behavior was given in 
(3.25 and 3.27)). . 

Let us denote the two terms on the right-hand side 
of (4.1) as follows 

Ma * ^AT-T’ ,meas ( 4 * 3a > h m A (fl)3£l < 4 * 3b ) 


Note that the derivatives have to be evaluated at a 

which Is unknown. Therefore, as usual, we use a 
Instead. This is based on the assumption that dfl Is 

small enough such that fl is close enough to fl. Define 


Gilfl) 


>A( fl )| 


(4.15) 


then using (3.1) we obtain 


—A 

01 

A 

02 

A — 

03 

A 

02 

A 

-01 

A J 

04 

A 

-03 

A 

-0 4 

A 1 

-01-1 


(4.16a) 


-0 2 

01 

-04 j 

A 

A 

A 

01 

02 

03 

A 

A 

A 

-04 

03 

-02-1 

(4.16b) 
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-A A 

-q 3 q 4 

d i ! ! 94 93 

1 

-O > 

rv> 
1 

r A A 

A J A A 

A ! 

G 3 " • ( >4 -93 

02 g 4 * 1 ’O 3 04 

91 

A * 

A j A A 

A 

LOj 02 

O 3 J l_02 "Oi 

94-1 

• 

(4.16c) 

(4. 16e) 

Further define 


hi > 

(4.17) 

Then (4.14) can 

be written as 


Y A - A(g)Yj - 

[ hi 1 h 2 1 h 3 1 h 4 ] |3qJ 1 
1 1 
|dq 2 | 

1 1 
Idq 3 l 
1 1 
|dq 4 | 

L J 

(4.18) 

Finally, define 

H < 

- C hi 1 h 2 1 h 3 1 h 4 ] 

(4.19) 

then (4.18) can 

be written as 



Y A - A(3)Yi - H q dg 

(4.20) 


How recall the definition of y A and ¥j as shown in 
(4.3). From these definitions it is obvious that we may 
substitute (4.11) and (4.20) into (4.1). When this is 
done and in view of (4.2), we obtain 


These two angles are A and B. Ui ,„ y 
quantities (tanA) m and {tanB) m , the 
measured by the sensor is computed as foil 


Using the measured 
unit vector 
ows 


lTtanA)"| 

I j 

-T’,meas - U + ( tan *)n + ( tanBJ " 1/2 | ( tanB) m | (4.25) 

L 1 J 

Let u m - (tanA) m and v m - (tanB) m then (4.25) becomes 


-T’ ,neas * + u m 2 + v m 2 )' 1/2 j v m 

IJ _ 


(4.26) 


Perturbation of (4,26) yields the following vector of 
errors for the measured sun vector. 


I du| 

dW T - - (1 + u„ 2 + v ra 2 )‘ ] /2 j dv j 


2. 


I 


u m du + u m v m dv ! 


' (1 + u m 2 + v m 2 )' 3/2 j Vm du + v m 2(iv 


u m du + v m dv -I 


(4.27) 


y ■ Hqdfl + M A y>QWj» + M A y » dWj » 
Note from (4.8) that 

(4.21) Let 

Q -d + u 2 

* vlr l/z 

9 * -[fix] 

(4.22) 

Wu - Q - 0 3 u 2 

W 12 * -Q 3u m v i 

therefore (4.21) can be written as 


W 21 ’ -Q 3u m v m 

w 22 - Q - q 3 > 

y - Hqdfl - M A ji [fixjy-p + M A y»dWj» 

(4.23) 

«31 - 

w 32 - -q 3 V[ 


The matrix M A t» Is not known to us; however, we do 
know Mat, It is easy to see that using the latter 
rather tnan the former does not affect the accuracy to 
any meaningful degree. For identical reasoning we use 
^T\mfias rat ^er Wj» * When these changes are made 
and' fhe order of the cross product is changed in 
(4.23) , we obtain 


Then (4,27) 


dy r - 


can be written as 

r«n w 12 -ifdu-| 
i hi 

I w 21 w 22 iLdv I 
I I 

L w 31 w 32-I 


I du | 
W s I ! 
Ldv | 


(4.28) 


(4.29) 


I 

I 1 
I 


H q da + M AT^T , ,meas x ^ + M AT d ^T’ 


(4.24) 


The measured quantities (tanA) m and (tanB) m can be 
written as m 


While (4.1) indicates how to generate the effective 
measurement * which updates the estimate, (4.24) 
indicates the linear relationship between y, the 
attitude errors, the misalignment errors of the sensor 
whose measurements are being used and dW T », the total 
error generated by the sensor. The derivation of (4.24) 
is the first stage in finding the measurement matrix, 
H, (defined in (2.8b)) for each of the sensors used 
onboard ERBS. In order to conclude the development 
which will yield those H matrices, we have to express 
dWy»1n terms of the error states of each sensor which 
constitute a part of £ shown in (2,9). This is done 
next. 

Fine Sun Sensor (F$$) Measurement Model 

The Sun sensor measures the tangents of the two 
angles of the vector from the spacecraft to the sun. 


(tanA) m - tanA + C A tanA + b A + n A (4.30a) 

(tan8) m - tanB + CgtanB + bg + n g (4.30b) 

where 

C A ,Cg * scale factor errors 

b A ,bR - biases modeled as Markov states In (3.25a) 
n A ,n B * white noise 

From (4.30) and the definition of u m and v m we realize 

+ h a f 


du - C A tanA + b A + n, (4.31a) 

dv - CgtanB 4 bg + n B (4.31b) 

When (4.31) Is substituted Into (4.29) the following is 
obtained 
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(4.38) 


dW- 


W|» 


w s l 


C A tanA 

l_CgtanB_ 


+ w c 


which can be written as 


I *>A i 

i i 


♦ W s l 


n A j 
_ n B-l 


(4.32) 


| | tanA 0 

dW T , - |W S | 

| | 0 tanB 


"II c a“ 

II 


- b B- 


+ W s | 




(4.33) 


This is then substituted into (4.24) resulting in 

1 1 0 ... 

.01 1 


1 • 1 H 1 0 . . . 

I_ 1 0 

• 0 j ^AT^-T’ ,meas x l j 


PtanA 

0 "II 1 0 

• 0 | 

M at W s | 

1 1 "atW s ! 0 

• 0 1 x 

| 0 

tanB || | 0 

• OJ 


i 

n A ! , 


♦ «AT«s! 

1 («. 

1 - 1 


l_ n B— I 


Equation (4,34) gives the measurement matrix, H, for 
the FSS which is used in computing the gain matrix and 
updating the covariance matrix. Since tanA and tanB are 
not available to us, we use, respectively, (tanA) m and 
(tanB)* Instead, Since the measured and the true 
quantities are close, this change practically 
introduces no error. 


IR Horizon Scanner Measurement Model 


The horizon scanner measures the roll and pitch of 
the spacecraft with respect to the geodetic coordinate 
system (GDS), i.e. it measures the direction of the 
nadir vector. The horizon scanner misalignment errors 
are assumed to be small with respect to roll and pitch 
errors, to be additive to roll and pitch and 
Indistinguishable from them. The unit vector in the 
direction of the nadir in the GDS is given as 

zJds - [0, o, 1] (4.35) 
In body coordinates this vector is given as 


| -cos (r)sin(p) 

I 

Zbody “ I sln(r) 

|_cos(r)cos(p) 


(4.36) 


where r is the roll angle and p is the pitch angle. 
As mentioned, this is the measured vector; that Is 


|-cos(r)s(n{p) | 

Mr.meas ’ | sin ( r > I < 4 - 37 > 

|_cos(r)cos(p)_| meas 

which is equal to the true vector plus error. The error 
vector is obtained by perturbing (4.36). The 
perturbation yields 


d Hr 


~sin(r)s1n(p)dr - cos(r)cos(p)dp 
cos(r)dr 

-sin(r)cos(p)dr - cos(r)sin(p)dp 


Let 


«h ■ 
then (4.38) 


| s1n(r)s1n(p) 

I 

I cos(r) 

I 

|_-s1n(r)cos(p) 
can be written as 


-cos(r) cos(p) 
0 

-cos(r) sin(p) 


dVjjt 


I J 

W h I dr | 

I J 
LdpJ 


(4.39) 


(4.40) 


We characterize the horizon scanner errors as bias 
(modeled as Markov process in (3.25b)) plus white 
noise; that is, 

dr* • d r + n^ r (4.41a) 

d P* - d p + n hp (4.41b) 

where d r and d p are the roll and pitch biases and n hr 
and n nD are the roll and pitch white measurement noise 
components. When (4.41) are substituted into (4.40), 
the following is obtained 


d *T\mea, ' Vh + w h0h < 4 - 42 ) 

where feu is as defined in (3.26) and nj * [ n h r * n bp^’ 
Since Is already in body coordinates, M AT given 

In (4.1) Ad in (4.24) is the Identity matrix. Since 
the horizon scanner was assumed not to have 
misalignment error the term containing misalignment 
angles in (4.24) is not needed. The model for the 
horizon scanner is given in (4.43) below. Again y is 
computed using (4.1). 

I" I 0 .... 0 | | 0 ... 0*| 

y - I H q I 0 .... 0 I w h I o ... 0 I* * W h n h (4.43) 

L q I o .... o I I 0 ... OJ 


Equation (4.43) yields the H matrix to be used with IR 
horizon scanner measurements . Similarly to the 
evaluation of the Sun sensor H matrix, we use the 
measured roll and pitch to evaluate W h in (4.43). 


Magnetometer Measurement Model 

The three magnetometers mounted orthogonally to one 
another measure the Earth’s magnetic field components 
along each of their axes. This arrangement of sensors 
Is identical to the three gyro arrangement which 
measure the spacecraft’s angular rate. The magnetometer 
error sources are also identical to the gyro error 
sources which are: scale factor errors, misalignments, 
bias (modeled as Markov process) and white measurement 
noise. Therefore the magnetometer errors can be 
represented by the same model as for the gyros. 
Therefore, In analogy to (3.13), we write the following 
expression for the errors Introduced by the 
magnetometers 
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| 3B x| | 5 mx 0 0_ |j B x"||5 Vy W|jVj 

j dB yj'[° V° || B yj + jV° VzjjBy H»Vy 

! dB zj |° 0 5 mz 1 1 B z I j 9 mzx Vy 0 j j ^z j j b mz I 


+ !v! 

(4.44) 

!».,] 

where B x , B y and B z are the magnetometer measurements and 

^my’ s mz^ 

(4.45a) 

* f 0 mxy» e mxz» e myx’ e myz» 0 mzx» e mzy^ 

(4.45b) 

* f^mx» d my» ^mz^ 

(4.45c) 

“ £ n mx» n my » n mz^ 

(4 . 45d) 


Substitute (4.46) into (4.51) and in view of (4.3b) 
also substitute (4.20) into (4.51). This results in 

X - MaT^T + M AT d S - A(fl)Yj + H q d 3 (4.59) 

Note that the first and third terms on the right-hand 
side of (4.59) cancel one another. Then when (4.50) is 
substituted into (4.59) we obtain the desired result 

1 - Hqdg + M AT B ’ x + + (4.60) 

or more expl icltly 

f I 0 0 | “| 

1 “ I « q | 0 0 | M at B’ | + M AT n m ( 4 . 61 ) 

|_ I o o I j 

Equation (4.61) yields the H matrix to be used with the 
magnetometer measurement updates. As for the previous 

sensors, we use the measured magnetometer data to 
evaluate B\ Finally, note again that the effective 
magnetometer measurements which have to be processed by 
the EKF are computed using (4.51) and not (4.1). 


5fl}» and fe- are the vectors of magnetometer scale 
factors, misalignments and biases, and n^ is the white 
measurement noise vector. The magnetometer output 
vector Bjnejs can be written as 

ST,meas " ^T + d S (4-46) 

where g T is the true magnetic field vector in the 
assumed magnetometer coordinates (note the difference 
between these sensors and the FSS and IR horizon 


scanner) . 

We 

can write 

(4.44) as follows 



|HB-| Tb x 

III 

0 

0 By B z 

0 0 

0 0 1 

0 

~0| 

1 


M 

I d By | • 1 0 
1 y \ 1 

B y 

0 0 0 

B x 8 Z 

0 0 0 

1 

0| 

j 

x + + 

1 1 
M 

|dB z | |0 

l:i l 

0 

B z 0 0 

0 0 

B x By 0 

0 

11 

J 


H 

where 








(4-47) 



* +T - 


®m» 




(4.48) 

Define the 

matrix B’ as 

follows 






1 

1 

CO* 

o 

o 

B z 0 

0 0 0 

1 

0 

0| 

j 


B’- 

10 

i 

By 0 0 

0 B x 

B z 0 0 

0 

1 

1 

0| 

(4.49) 


10 

L 

0 B z 0 

0 0 

0 B x By 

0 

0 

11 


We may write then (4.47) as 







dfi - B’a + + Dm (4.50) 


Measurement Error 


The main component in the magnetometer noise vector, 
Q-, is the quantization error. Its nature and 
characterization is explained as follows. The output of 
the magnetometers is received in the telemetry stream 
as 

N T - [N x , N y , N z ] (4.62) 

with H In counts. The lib component of N (i«x,y,z) is 
obtained from the actual measured components as [9] 

N i " INT ( B i,meas/ K i) (4-63) 

where INT means "the integer part of". Obviously, a 
certain part of the measured value is lost due to the 
INT operation: that is 

INT ( B i,meas/ K i) + "i * B i,meas/ K i «-6 4 ) 

The nature of the INT operation is such that can 

vary between 0 to 1. Moreover, the distribution of the 
chopped off value is uniform over the range 0 to 1. It 
is then easy to show that 

E(nj) - 0.5 (4.65a) Var{n i ) « 1/12 (4.65b) 

Substituting (4.63) into (4.64) yields 

N 1 ’ B 1,meas/ K 1 * "i (4.66) 

It is easy to see why in order to calculate the 
magnetometer readings on the ground the following 
computation is performed 


Note that dfi Is not identical to dft-p the way the 
latter Is defined In (4.4) since dfi contains the 
magnetometer misalignment errors whereas dy-r» does not. 
For this reason dfi cannot be substituted for dH r In 
(4.24). We have, then, to use the basic definition of y 
as applied to the magnetometer readings in order to 
formulate the linear relationship between y and the 
magnetometer errors which constitute For 

magnetometer measurements we formulate our effective 
measurement, y, as follows 


B i,comp * Kj[Nj + 0.5] (4.67) 

Substituting (4.66) Into (4.67) yields 

B 1,comp " K 1 t B 1 ,meas/ K 1 " n l + 0>5 ] (4.68) 

Define the measurement noise of magnetometer i as 

n m, 1 ‘ K l(°* 5 - "*) (4.69) 

then, in view of (4.65), 


I 

I 


“I 


* ■ M ATB T , me a S * MflHfi \ 


E{ n m, i) ‘ 0 (4.70a) Var(n mii ) - K?/12 (4.70b) 

where E denotes the expected value and Var denotes the 
variance. From (4.48) and (4.69) 
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B i,comp " B i,meas + n m,i 
and in vector form 

^comp " ^meas + 

From (4.46), (4.47), (4.61) and (4.72) it is obvious 
that the noise vector to be use in (4.61) is the vector 
r ^ as defined in (4.69) whose expected value and 
variance are defined in (4.70) 


(4.71) 

(4.72) 


V. QUATERNION NORMALIZATION 

The quaternion which represents attitude is a normal 
one. It was found [6] that forcing normalization on the 
estimated quaternion is advantageous since it speeds up 
convergence and eliminates the need for filter tuning. 
It was found in the present work too that normalization 
has these benefits. As shown in [6], normalization of 
the quaternion is equivalent to removing a portion of 
the estimate. This part that is removed must be 

accounted for in the next stage of the filtration. The 
handling of the normalization in this work is not 
identical to the one In [6] since here the covariance 
and state are propagated by solving their respective 
differential equations and not by using the state 
transition matrix as is the case in [6], There the part 
of the estimate which is removed by normalization is 
propagated using the state transition matrix and is 
considered at the next measurement update of the state 
estimate. Here though the normalization is done in 

between measurements. After the state estimate update 
by the horizon scanner measurement, the quaternion is 
normalized as follows 

3ip(+) * a j r (+)/ 1 3 1 r (+ ) I (5* 1 ) 

where the subscript IR denotes the fact that the 
quaternion estimate being dealt with is at the time 

point where the IR horizon scanner measurement are 
considered, (+) denotes the a*posteriori estimate and 
the superscript * denotes the resultant normal 
quaternion. It can be shown [6] that the normalization 

removes the following part from gj R (+) 

d<fy - a i r ( + ) flj R ( + x R (5*2) 

where dg IR - K r o y i r Is the estimate of dg which is 
computed using ihe kalman gain and the effective 
measurement of the scanner. Now when the FSS 

measurements are processed next, the estimate of the 
quaternion is updated as follows 

%$$(+) ■ d ^N + k fssI*f$s ' h fss^n 1 ( 5 - 3a ) 

ap$$(+) * Air( + ) + dflFSst*) ( 5 . 3 b) 

where flp$s(+) is the quaternion estimate after Its 
update by the FSS measurements. If no normalization Is 

performed, jjj R (+) ■ 3 ir(+K dg - 0 and (5.3) change 

accordingly. In any event, the quaternion, 2p$$( + )t is 
used as the a priori estimate of g for the magnetometer 
update, if available, or else Is propagated to the next 
time point. 


VI. COMPENSATION 

When propagating the state estimate and the 
covariance, we use the measured angular velocity. We 
know, however, that the propagated values are not 
accurate since the gyro outputs contain errors. As we 
estimate those errors, we can do better if we correct 
the gyro outputs for estimated errors. This operation 
is known as calibration. 


We also want to compensate the measurements obtained 
from the FSS, the IR horizon scanner, and the 
magnetometers which are all orientation measuring 
devices whose outputs are used to update the filter. 
The reason we want to compensate these sensors’ outputs 
Is different in nature than the reason for compensating 
the gyro outputs. Rewrite (4.1) and (2.11) 

* - "atKt\m.s - Mali! < 6 -U 

£ k (+) * ilcH + K k*k (6-2) 

The term KlYl is nothing but the estimate of x defined 
in (2.14) as 


X(t k ) - £{-) + X(t k ) (6.3) 

That is, in (2.11) we estimate the difference between 
the true value of X and its latest estimate, and add 
the estimate of the difference to the latest estimate 
of ^ to form its updated estimate. Now let us consider 
an error term in one of the sensor measurements, say a 
bias. This bias is a part of Wj> meas and thus, as 
Indicated in (6.1), bears its ’ signature on y. 
Consequently, if certain observability conditions are 
met, it is estimated and added to the state estimate as 
indicated in (6.3). If no compensation takes place, the 
next time the measurements of this sensor will be 
processed the bias will again be estimated and added to 
the previous estimate of this bias, thus creating a too 
large and hence wrong estimate. The correct way, then, 
to handle this case is to eliminate the estimate of the 
bias from Wy> meas - This way only residual bias which 
has not been ’estimated yet will be present in y as 
shown in (6.1). Only the estimate of this residual 
will, then, be added to the existing estimate of the 

bias, which is a part of £, yielding a correction to 
the previous estimate. This logic holds for the other 
error states too. The way we carried out the 
compensation is outlined in the ensuing. 


Gyro Compensation 




From (3.4) 

_w « w - 

dw 

(6.4a) 

therefore 





w - w - 

dw 

(6.4b) 

Rewrite (3.18) 




dw - 

[U|W| I] | 

'Sg + Ugl 

(6.5) 


1 

fig' 



l 

fcgj 



Since U and W are functions of w, and since the noise 
vector is of zero mean, a good estimate of dw is 
obtained from (6.5) as 


dw 


[U(fi)|W(w)|I] 



( 6 . 6 ) 


This estimate can then be used in (6.4b) to yield an 
estimate of y which is then used in the propagation 
algorithm instead of the raw gyro outputs. 

ESS Compensation 

Consider 
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Magnetometer Compensation 


* * M ATHT\meas 

- A(g)Vj 

(6.7) 

and recall (4.3) 



^A " ^AT^T’,meas (6.8a) 

¥a ‘ A (fl)*I 

(6.8b) 


As explained earlier, y is a linear function of £ which 
1$ the difference between & and We want * to go 
to zero when & approaches )(. Indeed, when & approaches 

<2 approaches g and Va computed in (6.8b) really 
yields the components of the unit vector in the Sun‘s 
direction when resolved in body coordinates. However, 

even when £ is equal to £, W-p meas will not be equal 
to Wj and therefore computed in (6.8a) will not be 
equaf to and thus y will not be zero. (This fact 
shouldn’t be confused with the fact that y goes to zero 
when the errors themselves go to zero). Only when 

Wji goes to W T as i goes to X, will y go to zero 

as It should. To achieve this use (4.4) to note that 

W T - M TT >(W T > irneas - dWf.) (6.9) 

As shown In (4.7), for small FSS misalignment angles 

Mj» j “1+0 (6.10) 

Transposing (6.10) we obtain 

M*j*j » ■ I - 0 (6.11) 

which when substituted into (6.9) yields 


As shown In (4.51) which is rewritten below, y is 
computed differently for the magnetometer measurements; 


namely 

* * H ATfiT,meas ' A (9)¥l 

(6.17) 

where 



and 

fiT.meas ’ Bj + di 

(6.18a) 


<12 - B’x + + Djj, 

(6.18b) 


x +T - [sj. aJ. tj] 

(6.18c) 


Following the rationale behind the FSS and IR horizon 
scanner and in view of (6.18), we compensate the 
magnetometer readings as follows. Compute 



where B’ is computed according to (4.49) using the 
uncompensated outputs of the magnetometers. Next 
compute the compensated magnetometer measurements 

- E T>me as * ( 6 - 19c > 

which are used to compute the effective measurement as 
follows 

1 m M^yly - A(g)Vj (6.19d) 

VII. IHE COMPLETE ERBS £KE ALGORITHM 


W T - (I - eHSr.rceas - d *r> < 6 - 12 > 

Therefore a reasonable estimate of Wy is 

0 T - (I - 9)(H T -, me as * d M < 6 ' 13 ) 

Replacinn Wy meas in (6.7) by Oy yields the desired 
result 


1 - M at (I - 0) (W T . >meas - dfl T .) - A(g)V, (6.14) 
where, in view of (4.32), dfif. is computed as follows 


| C A (tanA) m | I b A | 
dfi r - W s | a a I + W,l I 

S I _C b ( t a n B ) | 


(6.15) 


Note that initially when our estimate of £ is zero, 
(6.14) is reduced to (6.1). 


The models developed in the previous section were 
implemented into a program written in Fortran. The 
data used in the program is actual spacecraft data 
transmitted to Earth by ERBS. 

MEASUREMENT UPDATES 

The program is set up to compute an update initially 
and then propagate. The updates are performed for each 
sensor individually, the horizon scanner update is 
performed first, followed by the sun sensor, and the 
magnetometer. If any sensor data are not available the 
program bypasses that sensor and goes on to the next. 
In between each sensor update, the updated state and 
covariance are set to the a-priori values going into 
the next sensor update. If no sensor data are 
available, the a-posteriori state and covariance are 
set equal to the a-priori and are propagated to the 
next time point. 


IB Horizon Scanner Compensation 

The arguments made in outlining the FSS compensation 
are also valid for the horizon scanner only that here 
dMp is different and there are no sensor misalignments 
(they are considered to be a part of the bias errors). 
Following (4.42), we compute for the horizon scanner 

& v - 0 h £ h (6.16a) 

and A 

9 - 0 (6.16b) 

and substitute them into (6.14). The result is then the 
compensated effective measurement of the IR horizon 
sensor. 


Below is a summary of the algorithm and how it is 


applied to 

each update and to the propagation. 

IB Update 



Compute H: 

1 1 0 .. 

. 0 1 | 0 ... 0 

H - 

1 H 1 0 .. 

L 1 0 •• 

.01 w h 1 0 ... 0 

.01 1 0 ... 0_ 

Compute Wy’ 

,meas : 

|~cos(r)sin(p)~| 

j 


Bf’ ,meas 

- 1 sln(r) ( 

| ; 



|_cos(r)cos(p)_| meas 
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Compensate 


Magnetometer Update 


■ Hr.meas ' &h 

(7.3) 

Compute H: 


Compute residual and uncertainty of residual: 

* ■ - A(5)i I>1r 

(7.4) 

f 1 0 0 1 "1 

H - 1 H 1 0 0 | M at B’ | 

1 q 1 O 0 1 1 

(7.15) 

OJ 1r )j - {[HP(-)H t + R 1 r ])j ] / 2 

(7.5) 

Compute B r ,meas : 


where the subscript 1r denotes quantities pertinent to 
the horizon scanner, Vj is the spacecraft-to-earth 
unit vector obtained frdm the ephemerls and j denotes 
the jth element on the main diagonal of the residual 
covariance matrix. 

fB x "| 

Br.meas * jjyj 
Compensate fi T - ilneas : 

(7.16) 

State and Covariance Update 


P m ~! 


Compute K: 


’ BT’.meas ' | 

(7.17) 

K k - P k (-)H k T [H k P k (-)H k T 4 Rfc]- 1 

(7.6) 

IjJ 


Update £: 


Compute residual and uncertainty of residual: 


£ k ( + ) - £ k (-) 4 K k * k 

(7.7) 

* " At - A (3>^I,mag 

(7.18) 

Update P: 


(“mag)j * (CHP(-)H T + R^])^ 2 

(7.19) 


P k (+) - U-K k H k )P k (-)(I-K k H k ) T 4 K k R k K k T (7.8) 
SUB Sensor Update 


Compute H: 


| 0 .. 0 


~tanA 0 


H - | Hg | 0 . . 0 }fTieas x] M AT w s 

I ^ i n n 


II 


I 0 .. o I 


0 tanB_| | 


I 10 0 | 

I M at W $ | 0 0 | (7.9) 

I 5 I 0 0_| 


Compute M T ’ ) meas : 

Mr.meas ' 0 + ( tanA )m 2 + (tanB) m 2 ) ' ! / 2 
Compensate tfT’,neas : 


'(tanA)“! 
( tanB >fnj 
. 1 J 

... (7.10) 


where the subscript mag denotes quantities pertinent to 
the magnetometers, Vj is a unit vector in the 
direction of the magnetic field obtained from a 1980 
International Geomagnetic Reference Field model 
available In a Fortran subroutine and j denotes the jth 
element on the main diagonal of the residual covariance 
matrix. The state and covariance are then updated as in 
(7.6) through (7.8). 

STATE AND COVARIANCE PROPAGATION 

After all the sensor updates are performed, the 
state and covariance are propagated using a fourth 
order Runge-Kutta routine. The state and covariance 
are propagated ahead using the gyro data at the time of 
the update and the gyro data one second (nominally) 
ahead. Before propagating, though, the gyro data is 
compensated as follows. 


I ** 

s - [U|W| I] I L 


(7.20) 


St ■ (I - §)( Hr, me as - W s pA(t*" A >m"| ' w sjVj > 

... (7.11) 

Compute residual and uncertainty of residual: 

X - S T - A(fl)ii,fss (7 - 13) 

(u fss )j - ([HP(-)H T + Rf SS ]}j 1/Z ( 7 - 14 > 

where the subscript fss denotes quantities pertinent to 
the fine Sun sensor, 5ti f ss the spacecraft-to-Sun 
unit vector obtained frori a Solar-Lunar-Planetary file 
and j denotes the M element on the main diagonal of 
the residual covariance matrix. The state and 
covariance are then updated as in (7.6) through (7.8). 


State Propagation 


i(t) - £(i(t) ,t) (7.21) 

Covariance Propagation 


P(t) - F(i(t),t)P(t) + P(t)F T (£(t),t) 4 Q(t) (7.22) 


VIII. RESULTS 

Reference Solution 

The reference we used for comparison was the 
attitude solution obtained from the batch estimator 
used on the ground for operational attitude 
determination for ERBS. Table 8.1 shows the attitude 
solutions in the GDS and uncertainties for three 
different conditions (cases). Case 1 used sensor 
standard deviations of: FSS - 0.002 deg, IR » 0.2 deg, 
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MAG - 0.5 (unit vector). An orbit’s worth of data was 
used to compute the solution. Case 2 used the same 
amount of data as the first but the sensor standard 
deviations were: FSS - 0.1 deg, IR - 0.5 deg, MAG - 0.5 
(unit vector). Case 3 had the same sensor standard 
deviations as the second but only 30 seconds of data 
were used. Obviously, the reference s olution ii Ml 
unique although the real solution is. 


Table 1 Final Attitude Solutions 
IDd Uncertainties IdeaJ 


I 

| Yaw 
(Roll 
IPItch 


r 

Case 1 1 

Case 

i 

2 1 

. __ _ 1 . 

Case 

i 

3 1 

j 

Attitude 

Unc, i Attitude 

Unc. i Attitude 

iinsLJ 

1 




1 


-0.294 

0.01351 

-0.262 

1 

0.004 | 

-0.731 

1 

0.012| 

0.400 

0.01181 

0.421 

0.005) 

0.316 

0.070| 

0.650 

0.173 | 

0.420 

0.004 1 
1. 

0.384 

0.007| 

1 


Filler Solution 

Since ERBS is not inertially fixed, it is not very 
enlightening to see the variation in the quaternion. 
Therefore, in order to compare the filter solution to 
the batch solution, the estimated quaternion was 
converted to roll, pitch, and yaw in the GDS. We used 
the filter first to estimate attitude only. Figure 1 
shows the yaw solution in the GDS and is a typical 
example of the behavior when estimating attitude finly. 


YAW vs TIME 



Figure 2: YAW vs. TIME 

Full State 


YAW vs TIME 



Time (seconds) 


Figure 1: YAW vs. TIME 

No Calibration States 


The filter was then run with the full state starting 
at an a-priori attitude of zero degrees yaw, roll, and 
pitch. Figures 2 and 3 show the behavior of the yaw 
and pitch. Roll is similar to pitch. Figure 4 shows 
the estimation of the Z component of the gyro bias and 
Figure 5 shows an example of the residual behavior. 


PITCH vs TIME 



Figure 3: PITCH vs. TIME 

Full State 
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FSS Residual Y and Uncertainty , ''yro bias Z (deg/hr) 


GYRO BIAS Z vs TIME 



Figure 4: GYRO BIAS Z vs. TIME 

Full State 


FSS RESIDUAL AND UNCERTAINTY vs TIME 


0,06 


0.05 


0 . 0 * 


D.OJ 


0.02 


0 . 0 ! 


0 


- 0.01 


- Q 

1 

l 

□ FSS Residual 
Uncertainty 


. — 


0 20 *0 60 


Time (seconds) 


Table 8.2 gives the final attitude solutions and 
the various calibration states (after 60 sec of data) 
for 2 different sets of initial conditions (which 
actually differ only in the FSS and IR bias 
uncertainty) . 


Table 8.2: Attitude and Cal ibration States 
for Varying Initial Conditions 


Case 1 


Case 2 


| Yaw 

| -0.512 

i - 

i 

-0.676 


| Roll 

| 0.354 

i 

0.405 


| Pitch 

| 0.339 

i 

0.253 



i Value 

Unc. j 

Val ue 

Unc. 

IGYRO 

i 

i 

i 



[ SF 

| 0.194E-4 

0.01 | 

0 . 3 6 0 E - 5 

0.01 


| -0. 598E-4 

0.009981 

3.898E-5 

0.00998 


| -0.786E-4 

0.01 | 

-4 . 384 E - 5 

0.01 

|0 (deg) 

| 0.62E-4 

0.057 | 

3.2E-4 

0.057 


| 0.34E-4 

0.057 | 

0.50E-4 

0.057 


| 0.028E-4 

0.057 | 

0.013E-4 

0.057 


| 0.15E-4 

0.057 | 

0.046E-4 

0.057 


| -0.070E-4 

0.057 i 

0 . 052E-4 

0.057 


|-14. 6E-4 

0.057 | 

-10.2E-4 

0.057 

(Bias 

| -0.0211 

1.995 | 

-0.0984 

1.995 

| (Oeg/hr' 

| 0.0094 

1.995 | 

-0.0096 

1.995 


| 0.4381 

1.995 1 

0.3037 

1.995 

| ESS 

i 

i 

i 

i 



|0 

| 0.0648 

0.050 | 

0.0017 

0.056 

1 (Deg) 

! 0.0466 

0.054 | 

0.0114 

0.056 


| 0.0896 

0.057 | 

0.0140 

0.057 

| SF 

| -0.0235 

0.00607| 

-0.00557 

0.009 


| -0.0177 

0.00830| 

0.00044 

0.00954 

| Bias 

| 0.1732 

0.086 | 

0.077 

0.151 

j (Deg) 

I 0.0812 

0.058 | 

-0.0468 

0.0709 

1 IE 


1 



| Bias 

| 0.2115 

0.032 1 

0.0340 

0.0472 

1 (Deg) 

| -0.2377 

0.033 | 

-0.0840 

0.046 

IMS 

1 

i 

I 

i 



| Bias 

| -0.499 

0.699 | 

-0.597 

0.698 

1 (mg) 

| -1.479 

0.816 | 

-1.237 

0.816 


| 0.289 

0.911 | 

0.223 

0.911 


(other magnetometer states are negligible) 


Initial Uncertainties : 
Gyro 


Figure 5: FSS RESIDUALS Y vs. TIME 

Full State 


SF - 3*0.01 
0 - 6*0.057 deg 
bias - 2.0*3 deg/hr 
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FSS 

Q - 3*0.057 deg, 

SF - 2*0.01 

Bias « 2*0.1 in case 1, 
2*0.316 in case 2 

IR 

Bias - 2*0.1 in case 1, 
2*0.316 in case 2 

HAG 

Bias - 3*1 mGauss 


IX. CONCLUSIONS 

The ERBS EKF shows good, quick convergence 
properties when estimating only attitude. The filter 
is robust in that it can overcome initial attitude 
errors of up to 30 degrees (it may even go higher but 
30 degrees was the limit of our testing). When the 
remaining calibration states are added, and the sensor 
measurements are compensated for their calibration 
states the filter is not very robust. Starting the 
filter with a large initial attitude error would be 
outside of the linear region and the filter is not 
expected to give good behavior in those conditions. 

The attitude solutions estimated by the filter show 
some oscillation. Since the results presented are from 
real spacecraft data the filter would be expected to 
follow more closely the oscillations in the data 
whereas a batch solution would have these oscillations 
averaged out. 

We found, when estimating the entire state, that the 
results were dependent on the initial uncertainties due 
to a lack of observability. The batch solution which 
was used as our basis of comparison also was dependent 
on initial conditions. It could not be used as a true 
reference. 

The ability of the filter to quickly converge to an 
attitude solution from a large initial error 
demonstrates the feasibility of using an EKF for ground 
attitude processing in FDD, particularly in a real-time 
situation. Since all the states cannot be estimated 
simultaneously due to a lack of observability, more 
investigation into the a-priori uncertainties is 
necessary in order to achieve a desired accuracy in the 
final calibration states. 


X. FUTURE WORK 

In this work, the batch solution served only as a 
basis for comparison. It cannot be treated as a true 
reference. Simulated data will be used in the future, 
which will provide a true reference. From there 
further studies of the different calibration states and 
the ability of the filter to estimate them can be 
determined. At the time of this writing, the runs 
using simulated data were still being debugged. 

A further enhancement of the filter state will need 
to be made to include more than one sensor of a single 
type. Currently the filter only estimates calibration 
parameters for the sensors with coverage; no switching 
is done when the coverage changes over to another 
sensor of the same type. 


As mentioned previously, the ability to overcome 
large initial attitude errors makes the filter 
attractive for real-time operations. A real-time EKF 
will be developed which estimates attitude and possibly 
gyro calibration states only. 
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